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Context. With an average eccentricity of about 0.29, the eccentricity distribution of extrasolar planets is markedly different from the 
solar system. Among other scenarios considered, it has been proposed that eccentricity may grow through planet-disc interaction. 
Recently, it has been noticed that the thermodynamical state of the disc can significantly influence the migration properties of growing 
protoplanets. However, the evolution of planetary eccentricity in radiative discs has not been considered yet. 

Aims. In this paper we study the evolution of planets on eccentric orbits that are embedded in a three-dimensional viscous disc and 
analyse the disc's effect on the orbital evolution of the planet. 

Methods. We use the three-dimensional hydrodynamical code NIRVANA that includes full tensor viscosity and implicit radiation 
transport in the flux-limited diffusion approximation. The code uses the FARGO- algorithm to speed up the simulations. First we measure 
the torque and power exerted on the planet by the disc for fixed orbits, and then we let the planet start with initial eccentricity and 
evolve it in the disc. 

Results. For locally isothermal we confirm previous results and find eccentricity damping and inward migration for planetary cores. 
For low eccentricity (ej;2if/r) the damping is exponential, while for higher e it follows e oc e~ 2 . In the case of radiative discs, the 
planets experience an inward migration as long as its eccentricity lies above a certain threshold. After the damping of eccentricity 
cores with masses below 33 M Ea nh begin to migrate outward in radiative discs, while higher mass cores always migrate inward. For 
all planetary masses studied (up to 200 M Eanh ) we find eccentricity damping. 

Conclusions. In viscous discs the orbital eccentricity of embedded planets is damped during the evolution independent of the mass. 
Hence, planet-disc interaction does not seem to be a viable mechanism to explain the observed high eccentricity of exoplanets. 

Key words, accretion discs - planet formation - hydrodynamics - radiative transport - planet disc interactions - eccentricity 



1. Introduction 

One of the surprising orbital characteris tics of extrasolar plan - 
ets is their high mean eccentricity * 0.29 dUdrv & Santosll2007t) . 
Several explanations have been put forward to explain this dis- 
crepancy in comparison to the solar system. Planet-disc inter- 
actions are typically invoked to explain the planetary migration 
towards the central star that has occurred during the early for- 
mation phase. In addition to the change in semi-major axis, it is 
to be expected that the planet 's eccentricity will be influen ced 
through this process as well dGoldreich & Tremaind[l980h . It 
has then been suggested, by performing linear analysis, that 
the planetary eccentricity can be incre ased through the planet- 
disc interaction under some conditions (Goldreich & Sari 12003b 
ISari & Goldreichll2004t iMoorhead & Adamsll2008l) . Thev esti- 
mate that eccentric Lindblad resonances can cause eccentric- 
ity growth for gap-forming planets. Recently, a Kozai-type ef- 
fect between the disc and an inclined planet has been consid- 
ered dTerquem & Aimial2010l) . Numerical simulations, however, 
tend to show predominantly eccentricity damping for a vari - 
ety of masses dCresswell et afll2007l: IMoorhead & Ford1l2009]) . 
Additionally, the existence of resonant planetary systems with 
relatively low eccentricities (such as the system GJ 876) seems 
to indicate a damping action of t he disc on planetary eccentric- 
ity rather than an enha ncement dLee & Pealdl2002c iKlev et alJ 
l2005tlCrida et~al"1l2008h . 



On the other hand, very hi gh-mass planets can i nduce an 
eccentric instability in the disc dKlev & Dirksenll2006l) . In turn, 
the eccentric disc can possibly increase the plan etary eccentric- 
ity dPapaloizou et alj200lHD' Angelo et alj2006l) . However, this 
process can only explain the eccentricity of very massive (« 
5—10Mj U p) planets. Alternatively, planet-planet scattering seems 
to be a viable mechanism for increasing eccentricities through 
mutual gravitational interactions between the planets. The re- 
sulting eccen tricity distribution matches the observed one rea- 
sonably well dAdams & Laughlinl2003HJuric & Tremaind2008l : 
iFord & Rasiol2008l). Anoth er option is the fly-by of a nearby star 
dMalmberg & Daviesll2009l) . 

Planet-disc interactions have so far been studied mostly in 
the locally isothermal approach, where the temperature only 
depends on the distance from the central star. In this case 
for typical disc parameters, a negative torque is acting on the 
planet, and it migrates inward dTanaka et alJl2002l) . However, re- 
cently it has been shown that the torque acting on an embedded 
planet depends on the t hermodynamics of the dis c . Foll owing 
the pioneering work of IPaardekooper & Mellemal d2006l) . var- 
ious groups have now analysed the effect of the equation of 
state on the migration properties dBaruteau & Massed 120081: 
I Paardekooper & Papaloizoul l2008at Paardekooper & Mellemal 
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120081: iKlev&Cridal 120081) . Through full 3D radiative simula- 
tions of embedded planets, we have recently confirmed that 
including of radiation transport can produce a positive torque 
acting on low-m ass planets embedded in protoplanetary discs 
dKlev et al. 2009]), because through its action the required radial 
entropy gradient can be maintained in the disc. This results in 
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slowing down the inward migration, and under some conditions 
it may indeed be possible to reverse the inward migration pro- 
cess. 

The linear esti mates of the eccentricity evolut i on of 



embedded planets dArtvmowiczl 119931: IWard & Hahnl Il994t 
iTanaka & Wardl2004F concentrate on low eccentricities and pre- 



dict exponential decay on short timescales r e , 



where H/r is the aspect ratio of the disc and t, 



and 

T ecc the migration and eccentric i ty da mping timescale, re- 
spectively. IPapaloizou & Larwoodl (120001) have also considered 
larger values for e and they find an extended eccentricity 
damping timescale s uch th at de/dt oc e~ 2 if e > l.lH/r. 
ICresswell & Nelsonl d2006h have performed hydrodynamical 
simulations of embedded sma ll mass planets and find goo d 
agreement with the work by IPapaloizou" & Larwoodl J 2000). 
These 2D results have been confirmed bv lCresswell et alJ (120071) 
using fully 3D isothermal simulations. 

As mentioned above, the thermodynamics of protoplan- 
etary disc s is a crucial parameter for th e torque acting on 
the planet dPaardekooper & Mellemaf2006l) . Including radiation 
transport/cooling in a disc will give rise to positive torques act- 
ing on a planet embedded in such a disc, which indic ates out- 
ward migration. In our previous work dKlev et alJ2009t) we have 
shown that the inclusion of radiation transport/cooling in simula- 
tions with embedded low mass planets in three-dimensions (3D) 
can result in outward migration. So far these simulations have 
been limited to planets on fixed circular orbits. Now we extend 
this work and focus on the evolution of planets on eccentric or- 
bits for both, the isothermal and fully radiative regime. We focus 
first on low-mass planets and study the influence of the thermo- 
dynamics of the disc on eccentricity damping as well as on the 
evolution of the planet inside the disc. First we estimate the the- 
oretical migration and eccentricity damping rate for planets on 
fixed eccentric orbits. Secondly, we let the planets evolve in the 
disc and finally we investigate the influence of the planet mass 
on the change in eccentricity and semi-major axis. 



2. Physical modelling 

The protoplanetary disc is modelled as a three-dimensional 
(3D), non-self-gravitating gas whose motion is described by 
the Navier-Stokes equations. We treat the disc as a vis- 
cous medium, where the dissipative effects can then be de- 
scribed via the standard viscou s stress-tensor approach (eg. 
iMihalas & Weibel Miharaslll984l) . We also assume that the heat- 
ing of the disc occurs solely through internal viscous dissipa- 
tion and ignore the influence of additional energy sources (e.g. 
irradiation form the central star). This internally produced en- 
ergy is then radiatively diffused through the disc and eventually 
emitted from its surface. For th is process we use the flux-hmite d 
diffusion approximation (FLD, iLevermore & Pomraning|[T98ll) . 
which allows us to treat the transition form optically thick to thin 
regions approximately. The viscous forces used in our code are 
stated ex plicitly for the th ree-dimensional case in spherical coor- 
dinates in lTassoull (11978). We use a constant kinematic viscosity 
coefficient v with a dimensionless value of v = 1CT 5 (in code 
units, see below). This relates to the typically used ff-parameter 
through v = ac s H, where c s is the sound speed and H the ver- 
tical thickness of the disc. A more detailed prescription of the 
modelling and t he numerical me thodology is described in our 
previous paper dKlev et alJl2009h . We now extend the simula- 
tions, compared to our previous paper, by including planets on 
eccentric orbits with different masses. 



2. 1 . General Setup 

An important issue in modelling planetary dynamics in discs 
is the gravitational potential of the planet since this has to be 
artificially smoot hed to avoid singularities. We have shown in 
dKlev et al.ll2009l) that the physics of embedded planets can be 
described better by a cubic -potential rather than the often used 
e-potential. Hence, we use in this work the following form for 
the planetary potential throughout 
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Here nip is the planetary mass, d = |r— rp] denotes the distance of 
the disc element to the planet and r sm is the smoothing length of 
the potential measured in units of the Hill radius. The construc- 
tion of the planetary potential is in such a way that for distances 
larger than r sm the potential matches the correct l/r potential and 
is smoothed inside that radius (d < r sm ) by a cubic polynomial. 
The parameter r sm is equal to 0.5 in all our simulations, unless 
stated otherwise. 

The gravitational torques acting on the planet are calculated 
by integrating over the whole disc, where we apply a tapering 
function to exclude the inner parts of the Hill sphere of the 
planet. Specifically, we use the smooth (Fermi-type) function 
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which increases from at the planet location (d — 0) to 1 out- 
side d > Rh with a midpoint fi, = 1/2 at d = bRn, i.e. the 
quantity b denotes the torque-cutoff radius in units of the Hill 
radius. This torque-cutoff is necessary to avoid large, probably 
noisy contributions from the inner parts of the Roche lobe and 
to disregard mate rial that is gravitational bound to the planet 
dCrida et al.ll20Q9h . Here we assume (as in our previous paper) 
a transition radius of b — 0.8 Hill radii. 



2.2. Initial Setup 

The three-dimensional (r, 9, <p) computational domain consists 
of a complete annulus of the protoplanetary disc centred on 
the star, extending from r m ,„ = 0.4 to r max = 2.5 in units of 
fo = a j up - 5.2AU. In the vertical direction the annulus extends 
from the disc's midplane (at 9 = 90°) to 7° (or 6 = 83°) above 
the midplane for the simulations of planets with eccentric orbits 
in the midplane of the disc. Here 9 denotes the polar angle of 
our spherical polar coordinate system measured from the polar 
axis. The central star has one solar mass M* = M Q , and the total 
disc mass inside [r m j„,r max ] is Mdisc = 0.01M o . For the isother- 
mal simulations we assume an aspect ratio of H/r = 0.037 for 
the disc, in very close agreement with the fully radiative models 
of our previous studies. For the radiative models H/r is calcu- 
lated self-consistently from the equilibrium structure given by 
the viscous internal heating and radiative diffusion. The isother- 
mal models are initialised with constant temperatures on cylin- 
ders with a profile T(s) oc s with s = rsin#. This yields a 
constant ratio of the disc's vertical height H to the radius s. The 
initial vertical density stratification is approximately given by a 
Gaussian: 



p(r, 9) = p (r) exp 



(tt/2 - 9) 2 r 2 



2H 2 



(3) 
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Here, the density in the midplane is po(r) oc r~ ,5 which leads 
to a £(r) oc r~'^ 2 profile of the vertically integrated surface den- 
sity. In the radial and ^-direction we set the initial velocities to 
zero, while for the azimuthal component the initial velocity u$ is 
given by the equilibrium of gravity, centrifugal acceleration and 
the radial pressure gradient. This corresponds to the equilibrium 
configuration for a purely isothermal disc. 

For our fully radiative model we first run a 2D axisymmet- 
ric model (starting from the given isothermal equilibrium) to ob- 
tain a new self-consistent equilibrium where viscous heating bal- 
ances radiative transport/cooling from the surfaces. After reach- 
ing that equilibrium, we extend this model to a full 3D simula- 
tion, by expanding the grid into (^-direction. The resulting disc 
for this model has H/r « 0.037 so we choose that value for our 
isothermal runs. 



2.3. Numerical Setup 

Our coordinate system rotates at the initial orbital frequency 
of the planet (at r = ro). We use an equidistant grid in r, 0, <p 
with a resolution of {N r ,N g ,N,p) = (266,32,768) active cells 
for our simulations. At r m ,„ and r max we use damping bound- 
ary conditions for all three velocity components to minimise dis- 
turbances (wave reflections) at these boundaries. The velocities 
are relaxed towards their initial state on a timescale of approx- 
imately the local orbital period. The angular velocity is relaxed 
towards the Keplerian values, while the radial velocities at the in- 
ner and outer boundaries vanish. Reflecting boundary conditions 
are applied for the density and temperature in the radial direc- 
tions. We apply periodic boundary conditions for all variables in 
the azimuthal direction. In the vertical direction we set outflow 
boundary conditions for #„„■„ (the surface of the disc) and sym- 
metric boundary conditions at the disc's midplane (6 max = tt/2) . 
We use the finite volume code NIRVANA (IZiegler & Yorkdl 19971) 
with implicit radiative transport in the flux-limited di ffusion ap- 
proxim ation and the FARGO extension as described in lKlev et al.l 
(120091) . 



2.4. Simulation Setup 

In the first part of our model sequence we consider the orbital 
evolution of a planet with a fixed mass (20M Earth) on eccen- 
tric orbits using different initial eccentricities. For comparison 
we consider isothermal and fully radiative models. Using radia- 
tive discs here is a direct extension of a previous study under 
purely isothermal disc conditions using the same planet mass 
dCresswell et al. 2007). We distinguish two different approaches 
for these 20M£„ r ,/, models: first, a model sequence where the 
planet stays on a fixed eccentric orbit and secondly where the 
planet is free to move inside the computational domain under 
the action of the planet-disc gravitational forces. For the second 
models we insert the planet in the disc and let it move imme- 
diately, but using a time-dependent mass growth of the planet 
(through the planetary potential) until it reaches its destination 
mass. For the first set of models the IQMEanh planet is inserted 
as a whole in the disc at the start of the simulation. Initially the 
planet starts at a distance r - a j up = 5.2AU from the central star. 
For the fully radiative simulations we set the ambient tempera- 
ture to a fixed value of IQK at the disc surface (at #„„„), which 
ensures that all the internally generated energy is liberated freely 
at the disc's surface. This low temperature boundary condition 
works very well at optically thin boundarie s and does not ef- 
fect the inner parts of the optically thick disc dKlev & Linlll999t 



iKlev et al.ll2009h . In the second part of the project we consider 
sequences of models for a variety of planet masses. We note, 
that a lOMEanh planet has in our simulations using our standard 
resolution a Roche radius of about 3.3 grid cells. As we will see 
later in the results section, there is indication that for small mass 
planets and isothermal runs (using the cubic-potential) a higher 
resolution is required. 

3. Models with an embedded planet on fixed 
eccentric orbits 

In this section we consider planets remaining on fixed eccen- 
tric orbits embedded in either isothermal or fully radiative discs. 
From the disc forces acting on the planet we calculate its the- 
oretical migration rate and eccentricity change. Below we will 
compare this directly to moving planets in the isothermal and 
fully radiative regime. 



3. 1 . Torque and Power 

From the gravitational forces acting on the planet we can cal- 
culate the torque and energy loss (power) of the planet. These 
can be used to estimate the theoretical change of the eccen- 
tricity and the s emi-m ajor axis of the planet. Here we follow 
ICresswell etal] d2007h . 

The angular momentum L p of a planet on an eccentric orbit 
is given by 



L p = m p sjGM„a Vl - e 2 . 



(4) 



In our particular case, for non-inclined planets L p = L z . We can 
then obtain the rate of change of the semi-major axis and the 
eccentricity by 



1 a 

2 a 



disc 
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(5) 



where Tdisc is the total torque exerted by the disc onto the planet 



T disc = f (r p xF)\ z dV 

•J disc 



(6) 



Here r p denotes the radius vector from the star to the planet, F 
the (gravitational) force per unit volume between the planet and 
a disc element (at location r from the star), and dV the volume 
element. Equation (0) implies that a positive torque may also re- 
sult in eccentricity damping rather than outward migration. This 
is actually the observed result for our moving planets. 

The energy change per time (power) of the planet due to the 
work done by the gravitational forces of the disc is given by 



disc — I 1 p 
j disc 



r„ ■ F dV . 



(7) 



The energy of the planet depends only on the semi-major axis a 
of the planet and is given by 



1 GM.jn„ 

En = " 

p 2 a 



(8) 



We can now obtain for the energy loss and semi-major axis 
change 



E p a P di , 



\E P \ a \E p \ 



(9) 
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Fig. 1. Time dependence of the torque for different isothermal 
models with planets on fixed orbits. Top: Circular planets for dif- 
ferent values of the viscosity. The curves are ordered from high 
to small viscosity. For smaller viscosities the settling towards 
equilibrium takes longer. Bottom: Evolution for our standard 
viscosity (v = 1CT 5 ) for different eccentricities. Due to the large 
time variability of the torque for eccentric orbits these curves use 
sliding time averaged values with a window width of 1 period. 

From equation (O and equation (01 we can calculate the theoreti- 
cal change of the semi-major axis and the eccentricity for planets 
on a fixed eccentric orbit for the isothermal and fully radiative 
cases. We find for the theoretical change of the semi-major axis 



Pdis 



2a 



■ Pdis. 



a \Ep\ GM*m p 

and for the change of eccentricity 
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2 a 
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with our calculated a. 

Our simulations with planets on a fixed eccentric orbit fea- 
ture eccentricities ranging from <?o = 0.0125 to en. = 0.4. To in- 
vestigate a possible change in the orbital elements of a planet we 
first analyse the torques and power acting on the planet for fixed 
eccentric orbits. The time evolution of the torques for different 
models is displayed in Fig. Q] The top panel refers to planets 
on circular orbits for different viscosities, as quoted in the cap- 
tion. For the smallest viscosity (v = 1CT 7 ) the torque is unsat- 
urated and evolves through long period oscillations towards the 
equilibrium value. The timescale of the oscillation is compara- 
ble to the libration time of a particle near the edge of the horse- 
hoe region (see App. [A}. For larger viscosity the equilibration 
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Fig. 2. Torque (top) and power (bottom) acting on a planet on a 
fixed eccentric orbit in dependency of the eccentricity and ther- 
modynamics of the planet. The torque and the power have been 
averaged over 20 orbits, taken from t — 140 to t — 160 orbits for 
the isothermal (solid red line) and from t — 120 to t — 140 for 
the fully radiative simulation (dashed blue line). 



time becomes shorter as the the viscous diffusion time shortens. 
The resu lts are in very good agreement with existing 2D sim- 
ulatio ns dPaardekooper & Papaloizoull2009t iPaardekooper et all 
1201 Oh . and confirm clearly that viscosity is a necessary ingredi- 
ent for torque saturation. Additionally, this result indicates that 
the intrinsic (numerical) diffusivity is much smaller than that 
given by our standard physical viscosity. In the lower panel 
of Fig. Q] we display similar curves for eccentric planets and 
the standard viscosity v = 10~ 5 . Due to the strong variabil- 
ity of the torque on the o rbital timescale for eccentric orbits 
(see ICresswell et al.l (120071) and below) we display time aver- 
aged torques. 

In Fig. [2] we display the torque and power acting on such 
a planet for isothermal and fully radiative runs in dependency 
of the eccentricity of the planet. As only every 10th time step 
was plotted in our output file we averaged the torque over 20 or- 
bits to minimise the numerical fluctuations due to this procedure. 
For eccentricities smaller than 0.1 the planet experiences a neg- 
ative torque in the isothermal simulations, while for higher ec- 
centricities the planet feels a positive torque. The torque reaches 
a maximum at e — 0.175 and settles down to a nearly constant 
value for larger eccen tricities, which is in good agreement with 
ICresswell et al.l (120071) . The differences in the absolute values of 
their torque compared to ours have their origin in the aspect ratio 
H/ r of the disc. Our torque is generally higher, as a result of our 
lower H/r = 0.037 compared to their H/r — 0.05 disc. 
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Fig. 3. Theoretical rate of change of the semi-major axis (a /a) 
for eccentric planets on fixed orbits. For eccentricities smaller 
than about 0.025 we find a positive migration rate for the fully 
radiative disc, which indicates outward migration. This can ac- 
tually been observed for moving planets later in the text. 



The torque acting on the planet is in general higher for the 
fully radiative simulations than for the isothermal ones. This 
phenomenon was observed in various simulations in the past for 
planets staying on fixed circular orbits, so it is no surprise that the 
torque for the fully radiative simulation is higher for planets on 
eccentric orbits as well. For low eccentric planets the fully radia- 
tive simulation yield positive torques (and power) in agreement 
with our previous results for planets on purely circular orbits. On 
circular orbits, inclusion of radiation transport/cooling, gives rise 
to a positive torque, which implies outward migration in contrast 
to the isothermal simulations. The effect is caused by corotation 
region that gives rise to a positive contribution to the torque in 
the case of a positive entropy gradient. Interesting is the nar- 
row range in eccentricity of this outward migration. Already for 
e as 0.03 the direction of migration is directed inward again. This 
effect is caused by the spatially narrow region tha t contributes to 
the positive torque on the planet dKlev e t al. 2009). For larger ec- 
centricities the isothermal and radiative results match reasonably 
well, but the radiative results are always slightly larger than the 
isothermal ones. Different effects can contribute to this offset: 
the difference in sound speed, the spatially varying H/r for the 
radiative runs in contrast to the constant value for the isothermal 
runs, or a difference in the corotation torques. 

In the bottom diagram of Fig. [2] the displayed power of the 
planet follows the trend of the diagram of the torque acting on 
the planet with one big difference: the power of the planet is 
always negative for the isothermal simulations (implying inward 
migration), while it is positive for low eccentric planets in the 
fully radiative scheme. For small e the torque and power are very 
similar for both cases, since they must be identical for e — 0. 

The consequences for the inferred change in semi-major axis 
and eccentricity are displayed in Figs. [3] and [4] The theoret- 
ical migration rate (Eq. ITOb for planets on fixed eccentric or- 
bits (Fig. [3) reflects our assumptions. When a planet has a high 
initial eccentricity (e > 0.2) the migration rate is nearly con- 
stant and inward for the isothermal and fully radiative case, but 
as soon as the eccentricity gets damped to a value smaller than 
e = 0.2 the inward migration increases by a factor of 2 to 3. The 
fastest inward migration is seen for planets with an eccentricity 
of e as 0.125 for both thermodynamic cases. If the eccentricity 
evolves to lower values this rapid inward migration is slowed 
down in the isothermal case. In the fully radiative case this pro- 
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Fig. 4. Theoretical rate of change of the eccentricity (e/e) for 
eccentric planets on fixed orbits. 



cess is even stronger, so that planets with a very low eccentricity 
(e < 0.025) have a positive migration rate, indicating outward 
migration. The positive migration rate is a consequence of the 
positive torque and power acting on the planet. This confirms 
very well our previous work of low-mass plan ets on circular or - 
bits in fully radiative discs migrating outward dKlev et al.ll2.009h . 
The causes for this outward migration are the same as for circular 
orbits, as we will see later on. The migration rate in the isother- 
mal case is faster for low eccentric planets (with e « 0.125) and 
a little bit slower for high eccentric planets compared to the zero 
eccentricity case. 

The theoretical damping of eccentricity in the isothermal 
case (Fig. |4| indicates eccentricity damping for all values of e, 
in agreement with our results for moving planets (see Fig. Q~3] 
below). For low eccentricity planets (e < 0.10) the damping of 
eccentricity is much faster than for high eccentric planets. As 
soon as the eccentricity of high eccentric planets is damped to a 
low eccentricity, the damping of eccentricity becomes faster, but 
will then come to a constant value. In the fully radiative case, 
the damping of eccentricity is slower for low eccentric planets 
(e < 0.10) compared to the isothermal simulations but is nearly 
the same for high eccentric planets. This means the damping rate 
for eccentric planets is somewhat bigger in the fully radiative 
case compared to the isothermal case. For very small eccentric- 
ities we even get a positive value for the change of eccentricity 
(cut off in the Figure) in the fully radiative case. In the isother- 
mal case we note the opposite effect: for very small eccentricities 
the change of eccentricity is somewhat larger than for slightly 
higher eccentricities. This phenomenon might have its origin in 
a numeric feature: for small eccentricities the calculations of the 
migration rate and the change of eccentricity become very sensi- 
tive. Interestingly the slower change in eccentricity for the fully 
radiative simulations results in a lower final eccentricity. For ec- 
centricities smaller than about 0.025 we find a positive eccentric- 
ity change in the fully radiative case, which is actually not seen 
in our simulations of moving planets, and is a result of small 
numerical inaccuracy combined with the division by e in Eq.fTTI 



3.2. Torque analysis 

To understand the behaviour of the total torque in more detail we 
analyse now the space-time variation of the torque and power 
of the planet. For that purpose we introduce the radial torque 
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Fig. 5. Torque density acting on a planet on a fixed eccentric or- 
bit in dependency of the eccentricity of the planet embedded in 
an isothermal disc (top) and a fully radiative disc (bottom). The 
snapshots are taken at t = 150 Orbits for the isothermal simu- 
lations and at t — 100 Orbits for the fully radiative simulations. 
The planet is located at apoastron in all cases. 



Fig. 6. Torque density acting on a planet on a fixed eccentric 
orbit (e = 0.10) during the time of one planetary orbit in an 
isothermal disc (top) and in a fully radiative disc (bottom). One 
can clearly see that the torque changes according to the position 
of the planet to the central star, making it clear why the torques 
have to be averaged. 



density T(r), which is defined in such a way that the total torque 
T"" acting on the planet is given by 

T tot = J "™ T(r)dr. (12) 

The radial torque density has been a very useful tool to investi- 
gate the origin of the torques in our previous work on planets on 
fixed circular orbits. In Fig. |5]we display F(r) for a selection of 
our planets. All the snapshots were taken when the planet was 
located at apoastron. Note, that F(r) changes during the orbit, as 
the planet moves on an eccentric orbit. 

On can clearly see, that the major contribution to the torque 
originates at larger radii for larger eccentricities by construction, 
as the planets location is at apoastron and its distance from the 
central star is r — a + e, with a — 1 .0. For eccentricities smaller 
than e — 0.1 Fig. [5] suggests a negative total torque, while for 
larger eccentricities a positive torque is assumed, which can be 
clearly seen in Fig. [2] One might also argue, that the torque act- 
ing on the planet for eccentricities larger than e — 0.1 should 
be much higher than shown in Fig. [2] But keep in mind that the 
torque acting on the planet changes during one orbit as the planet 
moves on eccentric orbits, see Fig. [6] 

In Fig. [6] the planet has initially a distance of r = 1.1 to 
the central star (the planet is located at apoastron) and after half 
an orbit it is nearest to the central star (at r = 0.9, the planet 
is located at periastron) and moves then further away from the 



central star to r — 1.1 again. The motion of the planet in respect 
to the central star is the reason for the change in the F(r)-function 
with respect to time. 

In Fig. Q we display the surface density of planets moving 
on fixed eccentric orbits with e = 0.025, e = 0.05, e = 0.10 
and e = 0.20 in isothermal discs. These plots are taken at t — 
150 Orbits. For all surface densities displayed the planet is at 
apoastron, meaning the x value of the planet is -(a + e), with 
a - 1.0, while the y value of the planet is for all cases. 

Despite the fact that for eccentric orbits the density struc- 
ture and flow patterns appear and disappear periodically in phase 
with the orbit, one can see for the e = 0.025 case clearly two spi- 
ral waves exerted from the planet (one in the outer disc (r > r p ) 
and one in the inner disc (r < r p )). The spiral wave structure is 
comparable to the zero eccentric case. 

For higher eccentricities (e = 0.05 and e = 0.10) the outer 
spiral wave is more pronounced. At the time the snapshot was 
taken, the planet lies in apoastron where it moves more slowly 
through the gas, meaning that it is overtaken by disc matter on 
orbits that lie both inside and outside the planetary location. We 
also note a significant density enhancement close to the vicinity 
of the planet, which lies in front of the planet and will drag the 
planet forward, exerting a positive torque, see Fig. [5] The reason 
for this phenomenon lies in the flow lines, which are distorted by 
the planet's gravitational potential, and come to a focus in front 
of the planet. The corresponding azimuthally averaged density is 
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Fig. 7. Displayed are the surface density maps for planets on 
fixed eccentric orbits at t = 150 Orbits for isothermal simula- 
tions with H/r — 0.037. The plots feature eccentricities of e — 
0.025, e = 0.05, e = 0.10 and e = 0.20 from top to bottom. The 
planet is located at apoastron, meaning (x p ,y p ) = (-(a + e), 0), 
for each displayed eccentricity. 




Fig. 8. Azimuthally averaged surface density for planets on fixed 
eccentric orbits in the isothermal case. The displayed densities 
correspond to the surface maps displayed in the previous Fig. [7] 
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Fig. 9. Torque acting on the planet during the time of 2 planetary 
orbits for planets on fixed eccentric orbits in the isothermal case. 
The displayed torques correspond to the surface density maps 
displayed in Fig. [7] 



displayed in Fig. [8] Interestingly, for higher eccentricities there 
is no gap visible anymore in the averaged S-profile due to the 
complex structure of the spiral arms. These are actually the rea- 
son for the high torque density displayed in Fig. [6](for t = 150 
orbits). At periastron this effect is reversed, leading to a negative 
torque acting on the planet. For e = 0.20 we even see a stronger 
distortion in the density structure at apoastron. A more eccentric 
orbit of the planet will reduce the speed of it at apoastron more 
and leads thus to a more distorted density structure, compared to 
the zero eccentricity case. 

Comparing the torque density of the fully radiative simula- 
tion (see Fig. [5} with the isothermal torque density one notices 
only small differences for the different eccentricities. The torque 
acting on the planet is slightly larger for the fully radiative sim- 
ulation compared for the isothermal simulation. For smaller ec- 
centricities the reason for this phenomena lies in the spiral waves 
exerted from the planet. For very low eccentric planets, the argu- 
ments for a positive torque acting on the planet are the same as 
for planets mo ving on fixed circu lar orbits, which is explained in 
much detail in lKlev etall d2009h . 

In Fig. [6] we plot the change of torque density during one 
orbit for the fully radiative case for a planet on a fixed eccentric 
orbit with e = 0.1. The orbital evolution of the torque density 
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Fig. 10. Displayed are the surface density maps for planets on 
eccentric orbits on fixed orbits at t — 100 Orbits for fully radia- 
tive simulations. The plots feature eccentricities of e — 0.025, 
e = 0.05, e = 0.10 and e = 0.20 from top to bottom. The planet 
is located at apoastron, meaning (x p ,y p ) = (-(a + e), 0), for each 
displayed eccentricity. 



T(r) for the fully radiative simulation does not differ much from 
the isothermal one. One exception is the absolute value of the 
torque density at the location of the planet during the time of 
evolution in one orbit, which is higher for the isothermal case. 
The change of the torque density during one planetary orbit has 
the same reasons as for the isothermal case. In Fig. [9] we display 
the torque acting on planets on stationary eccentric orbits in the 
isothermal case. The eccentricities of the planet correspond to 
those shown in Fig.Qfor the surface density. One can clearly see 
that planets with a higher eccentricity (to about e * 0.10) show a 
higher amplitude in the torque acting on the planet, but for even 
higher eccentricities the amplitude is being reduced. One should 
also keep in mind that the averaged torque acting on the planet 
is negative for small eccentricities while it is positive for larger 
(e > 0.10) eccentricities. During periapses (f = 99.5 and 100.5) 
the planet also experiences a large energy loss (not displayed 
here), which follows in trend the same structure as the torque 
profile (e.g. the energy loss is greatest for e « 0.10). This leads 
to an enhanced inward migration for planets with an eccentricity 
around e ss 0.10, as can be seen in Fig. [3] for the calculated 
migration rate. 

The top figure in Fig.[l0]displays the surface density for an 
e = 0.025 planet on a fixed eccentric orbit in a fully radiative 
disc. The surface density distribution shown in this fi gure is very 
simila r to those found for fixed circular orbits (see iKlev et all 
120091) . For higher eccentricities the spiral wave structure can no 
longer be observed as clearly as for the low eccentricity case, no 
matter whether we are in the isothermal or fully radiative regime. 
The overall surface density structure for high eccentric planets 
for the fully radiative case matches nearly the corresponding 
structure in the isothermal case. As we calculate the torque act- 
ing on the planet due to the interaction with the disc material it 
is not surprising that the torques acting on the planet for high 
eccentric planets are very similar for the isothermal and fully 
radiative simulations. The corresponding azimuthally averaged 
density is displayed in Fig.QT] The profiles look very similar to 
those of the isothermal case, where for larger eccentricities the 
gap becomes invisible in the E-profile. 

In Fig.Q~2]we display the torque acting on planets on station- 
ary eccentric orbits in the fully radiative case. The eccentricities 
of the planet correspond to those shown in Fig. [10] for the sur- 
face density. One can clearly see that planets with a higher initial 
eccentricity (to about e « 0.10) show a higher amplitude in the 
torque acting on the planet, but for even higher eccentricities the 
amplitude is being reduced. One should also keep in mind that 
the averaged torque acting on the planet is positive for very small 
eccentricities (e < 0.025) and negative for 0.025 < e < 0.10 
while it is positive for larger (<? > 0.10) eccentricities. During 
periapses (t = 99.5 and 100.5) the planet also experiences a large 
energy loss (not displayed here), which follows in trend the same 
structure as the torque profile (e.g. the energy loss is greatest for 
e ~ 0.10). This leads to an enhanced inward migration for plan- 
ets with an eccentricity around e « 0.12, as can be seen in Fig. [3] 
for the calculated migration rate. 

4. Moving planets on initial eccentric orbits 

To study eccentricity damping of a planet embedded in a proto- 
planetary disc dynamically we now change our configuration in 
such a way, that the planet is able to evolve its orbit in time. In the 
beginning of the simulation we put the planet in the disc, and let 
the mass of the planet grow in such a way that the planet reaches 
its final mass of m p = 20M Earth at f = 10 orbits. This way we 
do not disturb the density distribution of the disc as much as by 
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Fig. 11. Azimuthally averaged surface density for planets on 
fixed eccentric orbits in the fully radiative case. The displayed 
densities correspond to the surface maps displayed in the previ- 
ous Fig.fTOl 
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Fig. 12. Torque acting on the planet during the time of 2 plane- 
tary orbits for planets on fixed eccentric orbits in the fully radia- 
tive case. The displayed torques correspond to the surface den- 
sity maps displayed in Fig.fTUl 



putting the planet with its full mass immediately into the unper- 
turbed disc. 



4. 1 . Isothermal discs 

The eccentricity evolution of various planets with individual 
starting eccentricity can be seen in Fig.Q~3] Planets with an initial 
eccentricity lower than about e « 0.10 experience an exponen- 
tial damping of eccentricity (as soon as the planet has reached 
its destined mass), while planets with larger initial eccentric- 
ity adopt initially a slower damping. As soon as planets with a 
higher initial eccentricity reach an eccentricity of about e * 0.10 
they undergo an exponential damping of eccentricity as well. For 
eccentricities 0.10 < e < 0.15 the damping of eccentricity is 
accelerated compared to the damping for higher eccentricities. 
This was expected as a result of our calculations of the change 
of eccentricities for planets on fixed eccentric orbits (see Fig. 
In the end all planets have settled to about the same value of ec- 
centricity (e * 0.02), independent of their starting eccentricity. 
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Fig. 13. Time evolution of eccentricity of planets embedded in 
isothermal discs with H/r = 0.037 and individual starting ec- 
centricity. The vertical dotted line at t — 10 indicates the time 
at which the planets have grown to their full mass. The black 
dashed lines indicate an exponential decay with a fitted r ecc = 29 
shifted to obtain a match for eo = 0.1, 0.2, 0.3 and 0.4. 



Us ing linear analysis for small eccentricities lTanaka & W ard 
(120041) find that the mean eccentricity change (averaged over one 
planetary orbit) is given by 



de/dt 0.780 



with the characteristic time 

*■ wave — 



^q 2 \ 1 / c s y* 

M, J \aQp 



(13) 



(14) 



where q denotes the mass ration between the planet and the star 
and the local surface density at the planetary orbit. For our 
2QM Earth planet at 5.2AU in the isothermal H/r = 0.037 disc we 
find the characteristic time to be t wave = 7.70 orbits, which gives 
an eccentricity damping time scale of about T ecc = t wave /0.78 = 
9.88 orbits. Apparently, this theoretically estimated decay time 
scale does not match the fitted value of T ecc = 29 as obtained 
from our numerical results. This exponential decay time is indi- 
cated by the black dashed lines in Fig.Q~3] As will be seen below, 
this strong difference occurs only for the isothermal case and is 
much reduced in the fully radiative models. In previous simu- 
lations on the evolution of eccentric planets in isothermal discs 
the agreement between th eory and numerics has been much bet- 
ter (ICresswell et al.ll2~007h . We attribute the difference to two ef- 
fects, a smaller H/r and the usage of the steeper more realistic 
cubic- potential instead of the more shallow e-potential that was 
used in lCresswell et al.l (120071) . Indeed, an additional run for the 
eo = 0.20 case with the e potential using r sm = 0.8 yields a fit- 
ting parameter of r ecc = 17. This is much closer to the theoretical 
computed value. 

But not only the eccentricity evolves in time on a moving 
planet, but also the semi-major axis of the planet as can be seen 
in Fig. [14] The reduction of the semi-major axis is in direct corre- 
lation with the damping of eccentricity. For planets with initially 
low eccentricity (lower than e = 0.12), we see a immediately re- 
duction of the semi-major axis on a fast rate. This rapid inward 
migration is then reduced to a very slow migration rate when the 
planets reach their final value of eccentricity (e » 0.02). Planets 
with higher initial eccentricity first migrate inward at a slower 
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Fig. 14. Time evolution of the semi-major axis of planets embed- 
ded in isothermal discs with H/r = 0.037 and individual starting 
eccentricity. 



rate as their low eccentric counterparts, but at the time their ec- 
centricity reaches the above mentioned point of e « 0.10 their 
migration rate changes and they undergo a rapid inward migra- 
tion, as their eccentricity damps exponentially. This rise in the 
migration rate was actually expected from the calculations of 
the theoretical migration rate for planets on fixed eccentric or- 
bits (see Fig. [3]). When the eccentricity damping then changes 
to a slower decay the migration rate becomes nearly linear and 
the planets move inward with a constant rate again. The planet 
starting with a zero eccentricity attains a low eccentricity dur- 
ing time and is also the planet migrating inward with the slowest 
speed. These results lead to the conclusion that planets on eccen- 
tric orbits do migrate inward at a slightly faster speed compared 
to their circular counterparts and that planets starting from initial 
higher eccentricity migrate inward faster than those with a lower 
initial eccentricity. 

But not only the initial eccentricity of the planet changes 
the evolution of the planet, but also the disc parameters. Our 
above shown simulations used a H/r — 0.037 in the isother- 
mal case. We now compare these results with isothermal simu- 
lations featuring H/r = 0.05. In Fig.[l5]we display the change 
of the semi-major axis and the eccentricity of lOMsarth planets 
in an isothermal H/r = 0.05 disc. The other disc parameters 
are the same as for the H/r = 0.037 planets. On the one hand 
the overall trend that the planet starting from a zero eccentricity 
orbit has the slowest inward migration is also true for the simu- 
lations with H/r = 0.05. On the other hand the inward migration 
in a H/r = 0.05 disc is faster for all initial eccentricities. In a 
H/r = 0.037 disc the planets with an initial eccentricity lower 
than eo = 0.20 end up with a very slow inward migration, while 
in the H/r = 0.05 case all planets end up with nearly the same 
migration rate, even the planets with an initial high eccentricity. 
These findings are in very good agreement with the results of 
ICresswell etaT] d2007l) . 

From our disc-data we estimate a theoretical exponential de- 
cay time for eccentricity damping timescale as T ecc = 32.95 
orbits. This value is larger than the previous one for the lower 
H/r = 0.037 case by a factor of (0.05/0.037) 4 . Numerically we 
obtain T ecc = 35, see the black dashed curves in Fig. [15] Hence, 
the numerical exponential decay in the simulations matches 
quite well to the theoretical damping rate. We attribute this better 
match to the smoothing effect of the higher pressure in the disc. 
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Fig. 15. Time evolution of the semi-major axis (top) and eccen- 
tricity (bottom) of a 20MEarth planet in an isothermal H/r — 0.05 
disc. The black dashed lines indicate an exponential decay of ec- 
centricity with T ecc - 35 for eo = 0.1, 0.2 and 0.3. 



4.2. Fully radiative discs 

If we now let the planet evolve its orbit with time for the fully ra- 
diative case, we obtain a very similar behaviour to the isothermal 
case for the evolution of eccentricity (see Fig. [Tot, but a quite dif- 
ferent behaviour concerning the semi-major axis evolution (see 
Fi g . ITTb . The damping of eccentricity for the fully radiative discs 
proceeds on a comparable timescale to the two previous results, 
and we find for the exponential behaviour T ecc = 25. 

The inclusion of radiation transport/cooling seems to have 
only a little effect on the damping of eccentricity. The damp- 
ing of eccentricity is somewhat slower in the fully radiative case 
compared to the isothermal simulations. For planets with ec- 
centricities lower than about e = 0.10 the damping of eccen- 
tricity follows approximately an exponential law. The eccentric- 
ity damping for high eccentric planets is first slower, until the 
eccentricity reaches a value of about e = 0.10 and is then in- 
creased to an exponential value (see the fit in Fig. IT6b. This rise 
in the speed of eccentricity damping is expected from the calcu- 
lation of the theoretical change of the migration rate for planets 
on fixed eccentric orbits (see Fig. [4|. The reduction rate of ec- 
centricity slows down when it reaches e * 0.05. In the end, as 
for the isothermal case, the eccentricity reaches the same value 
for all starting eccentricities, but the absolute value of eccentric- 
ity for the isothermal case (e « 0.02) is about a factor of five 
higher as for the fully radiative case (e as 0.004). The charac- 
teristic time for eccentricity damping in our fully radiative disc 
for the embedded 2QM Earth planet is t wave = 15.1 orbits, leading 
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Fig. 16. Time evolution of eccentricity for planets with individ- 
ual starting eccentricity for the fully radiative case. The black 
dashed lines indicate an exponential decay with T ecc = 25 for 
e = 0.1,0.2,0.3 and 0.4. 

to a eccentricity damping time scale of T ecc = 19.37 orbits, if 
we use H/r — 0.037 and an adiabatic sound speed. The black 
dashed lines in Fig. [16] indicate the exponential decay for our 
fitting, T ecc = 25. This is not exactly matching our theoretical 
results but the agreement is satisfactory. Note that as soon as 
the eccentricity of our simulations reaches e as 0.05 we observe 
a difference between our numerical results and the theoretical 
eccentricity damping. The damping in the simulations is much 
slower than the theoretical damping, as the planet stops inward 
migration at this point and stays on an orbit with nearly constant 
semi-major axis, which might slow down eccentricity damping. 

The density structure near the planet is smoothed for the fully 
radiative case. As the mass in the disc and near the planet is able 
to cool, it will give a smoother density profile. If the eccentricity 
is damped to a small enough value, the effects of radiation trans- 
port/cooling can set in, and will give rise to a positive torque, 
which results in outward migration (see plot of the semi-major 
axis Fig. [P71 i. As long as the eccentricity of the planet is higher 
than as 0.03 the planet will migrate inward, even for the fully ra- 
diative simulations, but with a rate smaller than the correspond- 
ing isothermal simulations. When the eccentricity reaches the 
critical value of e as 0.10 we observe a bump in the semi major 
axis of the planet. We expected this as the migration rate deter- 
mined by planets on fixed eccentric orbits (see Fig. [3} has its 
minimum there. 

Starting from a zero eccentricity planet, which will not gain 
much eccentricity during its evolution, on e can see direct out- 
ward migration, as predicted in our paper dKlev et al.ll2.009b for 
planets moving on circular orbits. Planets starting from a higher 
initial eccentricity will have to damp their eccentricity to a very 
low value to feel a positive torque acting on them and thus lead- 
ing them to outward migration. The eccentricity needed for out- 
ward migration is the same for all planets in our simulations 
(independent of the starting eccentricity) with a average value 
e as 0.03 for reversal. As soon as the e-damping is complete the 
planet experiences a positive torque and will migrate outward, 
even for the highest initial eccentricity in our simulations. 

The inclusion of radiation transfer/cooling in disc for embed- 
ded low mass planets on eccentric orbits will result in outward 
migration as soon as the eccentricity of the planet is damped to 
a value near the circular eccentricity. If the planet reaches low 
eccentric orbits which are nearly circular, the density structure 
will also become nearly like the one of a planet on a circular 
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Fig. 17. Time evolution of the semi-major axis of planets with 
individual starting eccentricity for the fully radiative case for 
20MEarth planets. 

orbit. The planet will then migrate outward as if it was on a cir- 
cular orbit. This outward migration will then be stopped as soon 
as the planet reaches regions in the disc where the density and 
temperature are too low to support outward migration. 

Looking at Fig. [17] it seems, that the outward migration is 
stopped for all planets at a certain critical radius in the disc. 
Even though at some point one might expect a termination of 
the outward migration depending on the local disc conditions, 
we note that in our case this feature appears to be a result of in- 
sufficient numerical resolution. For more information about this, 
see App.lAl 

5. Higher mass planets on eccentric orbits 

As was shown in many previous works the migration rate of a 
planet embedded in a protoplanetary disc does not only depend 
on the disc's str ucture and therm odynamics, but also on the plan- 
etary mass (e.g. lKlev et al.l2009h . In our previous work we found 
that planets up to about 33Mgarth experience a positive torque 
(which indicates outward migration), while higher mass planets 
experience a negative torque (indicating inward migration). It is 
now very interesting to investigate the evolution of planets with 
higher masses on eccentric orbits. In this chapter we focus on 
the eccentricity change and migration of planets ranging for 30 
to 200 earth masses in the isothermal and fully radiative regime. 
We focus here directly on migrating planets, and we do not pro- 
vide a torque analysis as we did for the IQMEanh planet. 

5. 1 . Isothermal discs 

In the isothermal regime using the cubic r sm = 0.5 potential may 
cause numerical problems (as described above) which become 
even more severe when the moving planets have a higher mass. 
The potential becomes just too deep for higher mass planets to 
consider our results as correct for the grid resolution we use for 
moving planets. A deeper planetary potential will give rise to a 
much higher density distribution near the planet which will be- 
come unrealistic in the isothermal case as the disc is not able to 
heat up upon compression. So we use in this section the common 
e-Potentia l with r sm = 0.8 f or planets with higher mass (as de- 
scribed in iKlev et al l (l2009h \ which will give us smoother and 
more realistic results in this case. Again, we let the planet reach 
its final mass during a time of 10 orbits. By this way the distur- 
bances in the disc are not as big as by inserting the planet with 
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Fig. 18. Time evolution of the semi-major axis for planets in an 
isothermal disc (H/r = 0.037) with 30, 50, 80, 100 and 200 Earth 
masses. In the top graph the planets have an initial eccentricity 
of eo = 0.10 and in the bottom plot the initial eccentricity is 
e = 0.40. 
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Fig. 19. Time evolution of the eccentricity for planets in an 
isothermal disc (H/r = 0.037) with 30, 50, 80, 100 and 200 Earth 
masses. In the top figure the planets have an initial eccentricity 
of eo = 0.10 and in the bottom figure the initial eccentricity is 
e = 0.40. 



its full mass at once. As the mass of the inserted planets becomes 
higher this feature becomes more and more important. 

In Fig. [18] we display the evolution of the semi-major axis 
for planets with 30, 50, 80, 100 and 200M Earth in the isother- 
mal case for eo = 0.10 and eo = 0.40. The high mass planets 
(M > 5QM Earth) migrate inward at the same rate, in contrast to 
the 3QM Earth planet, when starting with eo = 0.10. The planets 
seem all to migrate inward on a linear scale. However when the 
planets start with eo = 0.40 we observe a different picture; now 
only the 80, 100 and 2Q0M Earth planets migrate inward at the 
same speed and faster than the planets with lower mass. The mi- 
gration speed is faster in the beginning for the planets starting 
with eo = 0.10 compared to those starting with eo = 0.40, mean- 
ing that a high initial eccentricity tends to slow down the initial 
inward migration for planets with a higher mass (M > 30M£ flrt /,) 
which is in contrast to the results found for the 20Ms a rth planet. 
As the eccentricity is damped during time the migration speed 
for the planets starting with eo = 0.4 becomes faster than for 
the eo =0.1 planets. The observed bumps in the evolution of 
the semi-major axis occur for all planetary masses at the same 
time when we observe a change in the damping of eccentricity 
(displayed in the bottom figure of Fig. |T9Tl. 

In Fig. [19] we display the time evolution of eccentricity for 
planets starting with eo = 0.10 and eo = 0.40. The eccentricity 
drops immediately at the start of the simulation for the eo = 0.10 
case for all planetary masses, but the damping seems slowest 
for planets with higher mass, and in the end all high mass plan- 



ets (M > 50) end up with an eccentricity below e = 0.02. 
In the eo = 0.4 case the eccentricity damping sets in as soon 
as the planets have reached their final mass (after 10 orbits). 
The IQQMEarth planet experiences initially the fastest eccentric- 
ity damping until e « 0.22 and then the damping is slowed down 
to a nearly linear damping. The smaller planets follow in prin- 
cipal the same trend in eccentricity damping, only that the ini- 
tial damping is slower compared to the 200M£ fl ,,/, planet and 
terminates at a lower eccentricity. After the eccentricity reaches 
e < 0.17 the 100 and 200M Earth planet have a very similar sub- 
sequent eccentricity and semi-major evolution; at e < 0.15 the 
SOMEarth planet joins this evolution. 

The initial damping of eccentricity for planets with a high 
initial eccentricity (eo = 0.40) depends on the planetary mass, 
meaning the eccentricity for planets with higher mass is damped 
faster. This faster eccentricity damping is accompanied by a 
faster inward migration for higher mass planets. This trends 
seems to stop as soon as the planets have cleared their gap 
and the eccentricity damping and evolution of the semi-major 
axis is nearly the same for all the different high mass planets. 
Fig. [20] shows the azimuthally averaged surface density at the 
time of 200 planetary orbits for the two different initial starting 
eccentricities. The profiles after 200 orbits look very similar, the 
largest difference occurs for the 200MEarth model which displays 
clearly a wider gap and slower migration in the long run. For the 
models with 80 and 100 Earth masses it appears that the migra- 
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Fig. 20. Azimuthally averaged surface density after 200 orbits 
for isothermal models with different planet masses. Top: initial 
eo =0.1, and Bottom: initial eo = 0.4. 



tion is slightly faster for the high eccentric case. However, the 
averaged profile do not give a clear hint toward the cause. 
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Fig. 21. Time evolution of the semi-major axis (top) and eccen- 
tricity (bottom) of a 30M£ ar ,;, planet in a fully radiative disc. The 
eccentricity shrinks in the same fashion as for a lOMsarth planet, 
but about 50% faster. The semi-major axis increases also in the 
same trend, but the outward migration starts about 50% later than 
in the 20Me o , v /, case. The outward migration is then stopped at a 
radius slightly smaller than for a 20M£ £I(/ / 1 planet. 



5.2. Fully radiative discs 

In the fully radiative regime we can use our suggested cubic po- 
tential without problems for higher mass planets as the inclusion 
of radiation transport/cooling in a disc prevents a large density 
build-up near the planet, as the temperature near the planet rises 
and stops further mass accumulation. 

In Fig. [2]] we display the change of the semi-major axis and 
eccentricity over time of a 3QM Earth planet for different initial ec- 
centricities. One can clearly see that the outward migration starts 
when the eccentricity is damped to a small value as we already 
expected from our results of the IQMEanh planet. The damping 
of eccentricity is about 50% faster than for the IQMEanh planet. 
This speed up in the damping of eccentricity is due to the in- 
crease of the planet's mass. We also observe a change in the 
damping rate of the eccentricity as soon as the planet reaches 
e « 0.10, as for our previous simulations. For an initial low ec- 
centricity (eo = 0.05) we observe an earlier outward migration 
until the planets semi-major axis reaches the aforementioned 
barrier of outward migration in our discs. This barrier is depen- 
dent on the planets mass, as the final semi-major axis is slightly 
smaller for the 3QME ar th planet compared to the 2QME ar th planet 
(see App.lAt. 

In our previous work we obtained a limiting planet mass of 
about 33 M Earth for the occurrence of outward migration (see 
iKlev et al.ll2009l) . This implies that planets with a higher mass 



will not migrate outward but inward in a fully radiative disc. In 
Fig. [22] we display the evolution of the semi-major axis of plan- 
ets with 20, 30, 50, 80, 100 and 200 Earth masses for planets 
with an initial eccentricity of eo = 0.10 and eo = 0.40. In Fig. [23] 
we display the eccentricity for the same set of parameters. 

For the planets with 20 and 30 Earth masses we observe 
outward migration as soon as the eccentricity is damped to a 
small value (e < 0.02). As this damping is only achieved for the 
eo =0.1 simulations during the run-time of our simulations, we 
can not see outward migration in the eo = 0.40 figure, but we 
have seen it for the 2QM Earth planet in Fig. [171 

The planets with 50, 80, 100 and 200 Earth masses migrate 
always inward, independent of the initial eccentricity value. 
However the inward migration is much slower for the SOMEanh 
planet compared to the high mass counterparts. The 80, 100 and 
IQQMEarth planets migrate inward with the fastest rate, but the 
relative speed of inward migration for these three planets does 
not differ that much as it does for the 50 and 80M Earth- If high 
mass planets (M > 50M£ a ,-,/,) have an initial higher eccentric- 
ity they migrate inward a little bit slower than their low eccen- 
tric counterparts. As an eccentric orbit disrupts the typical spiral 
wave structure of a circular orbit the migration rate of a planet 
on an eccentric orbit is altered compared to the migration rate 
of a planet on a nearly circular orbit. Having reached a nearly 
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Fig. 22. Time evolution of the semi-major axis for planets in 
a fully radiative disc with 20, 30, 50, 80, 100 and 200 Earth 
masses. In the top figure the planets have an initial eccentricity 
of eo = 0.10 and in the bottom figure the initial eccentricity is 
e = 0.40. 
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Fig. 23. Time evolution of the eccentricity for planets in a fully 
radiative with 20, 30, 50, 80, 100 and 200 Earth masses. In the 
top figure the planets have an initial eccentricity of eo = 0.10 
and in the bottom figure the initial eccentricity is eo = 0.40. 



circular orbit the effects of radiation transport/cooling set in and 
the planet can (if its mass is low enough) migrate outward. 

In Fig. [23] we display the eccentricity evolution for planets 
starting with eo = 0.1 and eo = 0.4. The eccentricity damping 
for planets starting with eo = 0.10 seems to be independent of 
the planet's mass, if Mp > 25 MEarth - The damping sets in imme- 
diately after the planets reach their final mass and the eccentric- 
ity is damped to nearly the same value for all simulations. The 
planet's mass seems to have little effect on eccentricity damping 
for planets with a low initial eccentricity. The 100 and 200 A/ Earth 
planet show little fluctuations in the eccentricity when the eccen- 
tricity reaches zero. These fluctuations have their origin in the 
fact that this planet lies very close to the inner boundary of our 
simulation (we use reflective boundary conditions at the inner 
boundary) and thus interacts with it which give rise to the fluc- 
tuations in the eccentricity. These fluctuations are also enhanced 
by the planet's mass. 

On the other hand, if the planets have an initial high eccen- 
tricity (e.g. eo = 0.40), we observe a quite different damping 
rate of the eccentricity. The damping is faster for planets with 
higher mass, but the final eccentricity reached is the same for 
all planetary masses in our simulations. As soon as the 100 and 
IQQMEanh planet reach an eccentricity of e « 0.3 they loose half 
their eccentricity in a time of a few orbits, which also affects the 
migration rate of the planet, as we can observe a little bump at 
the same time. As in the isothermal case the initial damping rate 
of eccentricity is reduced for smaller mass planets. Interestingly, 



all planets experience a similar e-rate once the eccentricity has 
dropped below about 0.1 - 0.15. 

As for planets with 20M Earth an initial eccentricity influences 
the migration of planets with higher masses. In case the planet 
is small enough to undergo outward migration the effect of an 
initial eccentricity is a halting of outward migration to the point 
the eccentricity is damped to a very small value so that the effects 
of radiation transport/cooling can take effect as if the planet was 
moving on a circular orbit. For planets with a mass so high that 
they will not undergo outward migration an initial eccentricity 
can slow down the inward migration process by only a very small 
amount. 

5.3. The criterion for outward migration 

In the previous simulations we have seen that planets can migrate 
outward only for sufficiently small orbital eccentricities. The oc- 
currence of outward migration is linked to the detailed structure 
of the horse-shoe region since the to rques originate fr om a re- 
gion with a very small radial extent dKlev et al.ll2009t) . It is to 
be expected that eccentric orbits will destroy the detailed bal- 
ance, and this is what we see in our simulations. Nevertheless, it 
is interesting to estimate the value of eccentricity at which this 
reversal of migration can occur. For that purpose we have run ad- 
ditional series of test simulations in only two spatial dimensions 
but with identical physical disc parameter for various planetary 
masses. We measured the limiting value of the eccentricity with 
two alternative methods. In the first set we performed Simula- 
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Fig. 24. The critical eccentricity for reversal of migration ob- 
tained using 2D and 3D simulations for fully radiative discs. 
Only planets having eccentricities below this curve are prone to 
outward migration. 

tions with planets on fixed orbits for different eccentricities and 
masses. The point of migration reversal (equivalent to a sign- 
reversal of the power) has then been determined from the time 
averaged torque and power measured after 100 orbits. The time 
averaging has been done over 5 orbits. In the second alternative 
we followed the orbital evolution of the planet in the disc start- 
ing with an initial eccentricity of eo = 0.10. As demonstrated 
above all planets reduce this initial eccentricity and will migrate 
outward in a radiative disc after the eccentricity has dropped be- 
low a certain value. From the time evolution we determine this 
critical eccentricity. This second method is used for our full 3D 
disc as well. 

The results obtained using these procedures are displayed 
in Fig. [24] for the three sets of simulations. While the general 
trend is similar in all three series there are nevertheless differ- 
ences. The 2D results have an average value of e cr „ = 0.027 
with a drop for smaller mass planets. The values are systemati- 
cally larger than those of the 3D runs with e cr u = 0.015. This ef- 
fect may be caused by the slightly higher temperature in the disc 
([H/r]2D ~ 0.045) compared to the 3D runs or due to genuine 
flow differences in 2D and 3D geometry. Due to this larger H/r 
planets may experience o utward migration fo r higher masses 
as well in a 2D geometry (iKlev & Cridall2008l) . Within the 2D 
runs the results obtained for fixed stationary orbits yield slightly 
larger values for smaller masses and drop below the results 
for ev olving planets for masses above 25 ME ar m- In IKlev et al.l 
d2009l) we have observed that the location of the torque maxi- 
mum responsible for the outward migration lies slightly offset 
of the planet location at r ~ 0.984. Interestingly, for an eccen- 
tricity of e — 0.015 the apoastron lies close to this position, and 
we expect a destruction of this effect. So the planet can only 
migrate outward if its eccentricity is smaller than the distance to 
the offset in the torque distribution, which explains why the fully 
radiative effects do not turn the planet to outward migration if it 
has a high eccentricity. This eccentricity has to be damped below 
this value before the planet can undergo outward migration. The 
offset dista nce appears to be relatively insensitive to the mass of 
the planet dKlev et alJl2009t) and does no t scale directly with the 
planet's Hill radius. iBaruteau & Massetl d2008l) even suggest a 
torque maximum directly at the location of the planet. Our sim- 
ulations seem to indicate a dependence on the thermodynamic 
structure of the disc, such as radiative diffusion and temperature 
(entropy) gradients. This issue still requires resolution. 



6. Summary and Conclusions 

We have performed full 3D radiation hydrodynamical simula- 
tions of accretions discs with embedded planets of different 
masses on eccentric orbits. In a first sequence of simulations 
we have analysed in detail the dynamics of a planet with a 
given mass of 20M£ C „,/, for the isothermal as well as fully ra- 
diative case. In the isothermal situation we studied the cases 
H/r = 0.05 (a value often used in planet-disc simulations) 
and H/r = 0.037 (the value that matches the fully radiative 
case). In both cases we find similar behaviour for the eccentric- 
ity and_s£mi;major axis e volution, and results that match those 
of lCresswelleTaLl (120071) very well. Small eccentricities (with 
e'SlH/r) are damped exponentiall y with a time scale g iven ap- 
proximately by the linear results dWard & Hahnlfl994l) . Larger 
eccentriciti es are damped initially a c cordin g to e <x e~ 2 in agree- 
ment with iPapaloizou & Larwoodl (120001) and ICresswell et al.l 
(120071) . The final value of the eccentricity does not depend on 
the initial eccentricity of the planet. 

The planet migrates inward in the isothermal regime. Low 
mass planets (e.g. 20MEarth) on eccentric orbits with large ec- 
centricity (e > 0.20) have a slower migration rate as their low 
eccentric counterparts in the isothermal regime. As soon as the 
damping of eccentricity proceeds to smaller values the migration 
rate is pumped up as if the planet had started with a low eccen- 
tricity. The maximum inward migration rate occurs at e as 2H/ r. 
In the fully radiative regime high eccentric planets (e > 0.20) 
with lOMEarth migrate inward on a rate comparable with the 
isothermal regime. The corresponding eccentricity damping rate 
for the fully radiative scheme is about the same as for the isother- 
mal simulations, taking into account the different sound speeds. 

But as soon as the eccentricity becomes small enough the 
planets undergo a change in the direction of migration. The in- 
clusion of radiation transport/cooling in discs with embedded 
low mass planets will give rise to a change in the direction of 
migration for planets whose initial eccentricity has been damped 
to a nearly circular orbit. The maximum eccentricity a planet 
can have to still undergo outward migration seems to be deter- 
mined by the torque maximum in our F(r) function. This torque 
maximum has a slight offse t comp ared to the planets location as 
demonstrated in IKlev et all d2009l) . corresponding to a limiting 
eccentricity of about 0.015-0.025. If the eccentricity of the planet 
is larger than this value it will migrate inward, while it will mi- 
grate outward if its eccentricity is smaller (see Fig.[24t. For very 
small planet masses the maximum eccentricity is reduced. For 
planets on nearly circular orbits the effects of radiation keep the 
positive torques acting on the planet unsaturated, which implies 
continuous outward migration. Moving planets experience this 
result directly, and do indeed migrate outward in the disc in con- 
trast to planets in the isothermal regime. 

Eccentric planets with higher mass are slowed down in their 
migration at the beginning if they have a high initial eccentric- 
ity (eo ^ 0.20) in the isothermal as well as in the fully radia- 
tive scheme. If e > 0.02 all planets move inward independent 
of their mass, even those embedded in a fully radiative disc. 
When the eccentricity is damped further the high mass planets 
(M > 50M Earth) still move inward for both regimes, as they open 
a gap in the disc and migrate as Type EL The eccentricity damped 
30MEarth planet moves outward in the fully radiative sc heme as 
we expected from our previous results dKlev et al.ll2009t) . 

Independent of the discs thermodynamics and planet mass, 
we find that an embedded planet with a given initial eccentricity 
will lose this eccentricity in time. The rate of the eccentricity loss 
depends on the value of the initial eccentricity but is much faster 
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than the migration time. Hence, according to our results planet- 
disc interaction cannot be the cause of the observed high mean 
eccentricity of extrasolar planets. This finding is supported by 
the fact that the existence of planetary systems in mean-motion 
resona nce with small libra t ion angles require damping of eccen- 
tricity dLee & Pealell2002t iKlev et alj|2005l) . A solution to this 
problem might be planet-planet interaction, which we have not 
considered here. 

While performing our studies we noticed that numerical res- 
olution is a serious issue in these type of simulations. As shown 
in the appendix, only for very high numerical resolution or in 
an adapted rotating coordinate system which rotates with the 
present location of the planet, can we observe the outward mi- 
gration. Finally, as the origin of this outward migration for plan- 
ets below M p « 33ME art h is created by a delicate balance of 
torques which is destroyed by even a very small eccentricity of 
the planet, the question arises how this effect can persist under 
realistic conditions. It remains to be studied what effect the tur- 
bulent motions within the disc have on the corotation torque of 
the planet. 
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Appendix A: Numerical features for moving Planets 

In Fig. [17] we have noticed that the outward migration of the 
moving 20 M Earth planet terminates at a well defined radius inde- 
pendent of its initial eccentricity. To determine theoretically the 
extent of outward migration from our disc properties it seems 
useful to compare different timescales: the libration time and 
the radiation time scale of the disc. The necessary unsaturated 
torques needed for sustained outward migra tion require approx- 
imatel y equal libration and cooling time (seelBaruteau & Massel 
(120081) : IPaardekooper & Papaloizoul (l2008bh N ). For the latter we 
use in our case the radiative diffusion time scale. We define 
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The libration time is (as in iBaruteau & Massel d2008l) ) t/,;, = 
%jirpl{3>Q.KX s ) with the half-width of the horseshoe-orbit x s = 
\.\6rp ^jql(Hlr). To compute the timescales we use the den- 
sity and temperature at the midplane of the disc at the begin- 
ning of the simulations. The timescales are displayed in Fig. IA.ll 
The timescales seem to be comparable to about r — 1 .4, so the 
planet should be able to migrate outward at least to this point and 
should not stop at r — 1.23 as observed in Fig. [171 In the plot we 
also show additionally the viscous timescale t v ; sc - s 2 /v which 
should be comparable to the radiative time for accretion discs in 
equilibrium. Apparently, this relation is well fulfilled. 

The migration of the planet inside the disc is cal culated via 
the to rque acting on the planet. In our previous work (IKlev et alJ 
120091) we found that the calculated torque acting on the planet 
appears to be converged already for our standard resolution if the 
simulation is performed in a rotating coordinate system where 




Fig. A.l. Timescales in dependence of the distance from the cen- 
tral star. To compute the timescales we used the density and tem- 
perature of the midplane at the beginning of the fully radiative 
simulations. 



the planet is fixed for circular orbits. If the planet is allowed to 
move freely inside the disc, for example due to the disc's grav- 
itational force, this is no longer possible, as the planet changes 
its semi-major axis during time. If the planet moves to an orbit 
with larger radius the rotating frame will not be able to keep the 
planet at the same point in the grid, as the rotation frequency of 
the frame is not linked to the planet anymore. 

To investigate this, we performed a series of simulations with 
a lOMsarth planet on a fixed circular orbits placed at various dis- 
tances from the star. In each of the runs the coordinate system 
rotated with the orbital period of the planet such that the planet 
did not move through the grid. In Fig. IA.2l we display the torques 
acting on this planet in the fully radiative scheme. The planets 
are placed in a disc, corresponding to our standard model, at dif- 
ferent planetary radii ranging from r p = 0.8rj up to r p = 2.0rj up . 
As the rotation frequency of the planet matches the rotation fre- 
quency of the coordinate system, the planet remains at a fixed 
location in the grid during the whole simulations. Clearly, the 
torque acting on these planets is positive for planets with r p up 
to r p - \.9rj up , indicating that a planet at these radii should still 
migrate outward. The positive torques are contrary to the results 
of the moving planets in Fig. [TTJ where the planets stop their 
outward migration already at about r — 1.23. 

The torques acting on theses planets indicate, that a planet in 
the disc should migrate outward until at least r = 1.9. Or possi- 
bly even further, if not stopped by our finite computational do- 
main where numerical effects from the outer boundary may dis- 
turb its way. To test this hypothesis we calculated the evolution 
of a 2QM Earth planet starting at r = 1.5 with a rotation frequency 
of the grid matching the rotation frequency of the planet at this 
distance (r = 1.5). Indeed, as Fig. IA.3I shows, the planet now 
migrates outward as the torques presented in Fig. lA.2l suggested. 
Note that the mass of the planet has been increased gradually 
within the first 10 orbits. Again, this result is in contradiction to 
the stopped outward migration in Fig. [TTJ 

May this effect be caused by the difference in the rotation fre- 
quency of the planet and the numerical grid? To answer this, we 
first place a 2QM Earth planet at r p = 1 .5 with a rotation frequency 
corresponding to a r = 1 .0 and let it evolve in the disc. The evo- 
lution of the semi-major axis is displayed in Fig. IA.4l Secondly, 
we started an identical planet at r p = 1.0 with the same rota- 
tion frequency of the numerical grid. In the end of the evolution 
both planets comes to a halt at the same radius. Obviously, the 
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Fig.A.2. Torque for planets on fixed circular orbits with dis- 
tances r p = Q.%rj up to r p = 2.0rj up . In each simulation the coor- 
dinate system is rotated with the planet. 
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Fig. A.4. Semi-major axis of planets starting at r — 1.0 and r 
1 .5 with a rotation frequency of the numerical grid matching the 
Keplerian rotation at r — 1 .0 for both cases. 
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Fig.A.5. Torque acting on a IQMEanh planet at r p - 1.5. Red 
solid line: The rotation frequency of the grid matching the 
Keplerian rotation frequency at r — 1 .0, and the grid resolution 
has been doubled compared to our standard simulations. Blue 
dotted line: Standard resolution with grid velocity equal to the 
planet. 
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Fig.A.3. Semi-major axis of a lOMsarth planet starting at r p = 
1.5 with a grid rotating with the initial angular speed of the 
planet. 



chosen rotation frequency of the numerical grid has an influence 
of the migration of the planet inside the disc. The reason lies 
in the increased numerical diffusion that occurs if matter moves 
with a fast speed through the disc. The implemented FARGO al- 
gorithm helps to solve this problem but cannot fully eliminate 
it. The problem is enhanced in our situation because it is the de- 
tailed fine structure in the flow near the planet that determines 
the outcome of migration. 

Being caused by numerical diffusion, the effect may be (par- 
tially) cured by increasing the resolution of the grid. We now 



double the grid size to 532 x 64 x 1532 active cells in r, 6, (p 
direction and let a planet fixed at r p = 1 .5 evolve with a rota- 
tion frequency corresponding to r = 1 .0. The torque acting on 
the planet is displayed in Fig. IA.5l and is clearly positive, indi- 
cating outward migration. Over-plotted is the torque of a planet 
at r p = 1 .5 with matching rotation frequency and our standard 
resolution (the same as in Fig. IA.2I >. 

Hence, it seems that the rotation frequency of the numerical 
grid influences the migration of the planet in the disc, if the nu- 
merical resolution of the grid is too small and the planet moved 
to a radius where its rotation frequency differs by more than 25% 
from rotation frequency of the grid. In that way, our results de- 
termined in the main article are correct, as the planets have mi- 
grated inside the disc only by a little bit before their eccentricity 
is damped and they start their outward migration. The obtained 
limit of r — 1 .23 for the outward migration seems to be a numeri- 
cal artefact, however. To obtain an accurate migration for planets 
under these conditions it seems best to perform the simulations 
in a coordinate system that rotates always with the speed of the 
planet. For multiple planetary systems this is not possible and a 
higher resolution is required. 
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